#include<iostream>
using namespace std;
#include<fstream>
#include<cmath>
#include<ctime>

#include"library.h"

extern "C" void dsyev_(char *a,char *b,int *n,double *A,int *nn,double *e,double *work,int *lwork,int *info);

int main(){

	int i,j,k,l, nkeep, n;
	double y;
	double *A, *a, *z, *q0;
	cout<<" matrix dimension = "; cin>>n;
	cout<<" nkeep = "; cin>>nkeep;

	srand(time(NULL));
	A = new double [n*n];
	for(i=0;i<n;i++){
		for(j=i;j<n;j++){
			A[i*n+j] = (double)rand()/RAND_MAX;
			A[j*n+i] = A[i*n+j];
		}
	}

	a = new double [n*n];
	for(i=0;i<n;i++){
		for(j=0;j<n;j++){
			a[i*n+j] = A[i*n+j];
		}
	}
	z = new double [n]; for(i=0;i<n;i++)z[i]=0;

	q0 = new double [n];
	for(i=0;i<n;i++) q0[i] = (double)rand()/RAND_MAX;
	y = sqrt( InnerProduct(n, q0, q0) );
	for(i=0;i<n;i++) q0[i] /= y;

	clock_t start = clock();
	Lanczos(n, A, nkeep, z, 1E-6);
	clock_t end = clock();
	cout<<"Lanczos takes "<<(end-start)/CLOCKS_PER_SEC <<" s.\n";

	char alpha='v', beta='L';
	int lwork = 3 * n, info=0;
        double *e = new double [ n ];
        double *work = new double [ 3*n ];

	start = clock();
	//dsyev_(&alpha, &beta, &n, a, &n, e, work, &lwork, &info);
	end = clock();
	cout<<"dsyev takes "<<(end-start)/CLOCKS_PER_SEC <<" s.\n";
	if(info<0){
		cout<<"The "<<abs(info)<<"th argument had an illegal value."<<endl;
		exit(1);
	}
	else if(info>0){
		cout<<"The algorithm failed to converge: "<<abs(info)<<" off-diagonal elements of an intermediate tridiagonal form did not converge to zero.\n";
		exit(1);
	}
	//cout<<"lapack, eig = "; for(i=0;i<n;i++)cout<<e[i]<<"\t"; cout<<endl;

	// for test only
	for(i=0;i<nkeep;i++){
		cout<<"A v / lambda:";
		for(j=0;j<n;j++){
			y=0;
			for(k=0;k<n;k++){
				y += a[j*n+k] * A[i*n+k];
			}
			if( fabs(y/z[i]) > 1E-4 )cout<<y/z[i]<<",";
		}
		cout<<endl;
		cout<<"v:"; 
		for(j=0;j<n;j++){
			if( fabs(A[i*n+j]) > 1E-4 )cout<<A[i*n+j]<<",";
		}       
		cout<<endl;
	}

	delete [] q0; delete [] z; delete [] a;
	delete [] A;

	return 0;
}
